* Creates dataset and then produces outputs for Table 4 (randomisation checks)

clear all

use "$savedata/masterdata.dta", clear

merge m:1 extract using "$savedata/past_op0517.dta"
keep if _merge==3
drop _merge

merge m:1 extract using "$savedata/pastae0717.dta"
keep if _merge==3
drop _merge

merge m:1 gpprac finyear using "$savedata/gp_pastadmits0517.dta"
drop if _merge==2
drop _merge

* Need to create "last year" variables for use of ED (A&E), OP, and cardio related visits - OP sample from 2007 onwards (eventually use 2009)
foreach var in num_cardio_op num_visit_op any_visit_op any_cardio_op{
	
gen prevyear_`var' = `var'2005 if finyear==2006

forval year=2007(1)2017{
	
	local past = `year'-1
	replace prevyear_`var' = `var'`past' if finyear==`year'
}

}

** ED sample only for 2009 onwards
foreach var in any_amb any_cardio_ae any_visit_ae num_amb_ae num_cardio_ae num_visit_ae{
	
gen prevyear_`var' = `var'2007 if finyear==2008

forval year=2009(1)2017{
	
	local past = `year'-1
	replace prevyear_`var' = `var'`past' if finyear==`year'
}

}

drop any_amb* any_visit* any_cardio* num_amb* num_visit* num_cardio*

keep if sample25==1
gen vol = vol25


merge m:1 lsoa using "$inputs/lsoa_areatypes.dta"
drop if _merge==2
drop _merge

* Do analysis from here

preserve

reghdfe survive365 c.prevyear_cost, absorb(i.sex##i.derv_age i.black i.mixed i.chinese i.asian i.race_miss i.ynch* i.prevyear_stroke i.di1 i.di2 i.di3 i.di4 i.di5 i.shock i.arythmia i.arthero i.arrest i.dow##i.admidate_mont##i.finyear i.hyid, savefe) keepsingleton
predict residual, residuals

xtreg residual c.std_ami3 c.std_nonami3, fe i(doctor_id)
predict docfe, u
gen sigma_u = e(sigma_u)
gen sigma_e = e(sigma_e)

gen signal_2step = ((sigma_u^2) / ((sigma_u^2) + ((sigma_e^2)/vol)))
gen adj_docfe = docfe*signal_2step

egen std_adj_docfe = std(adj_docfe)

foreach var in lsoa_imdscore10 lsoa_wa_benefits prevyear_admit_em prevyear_admit_el prevyear_los lsoa_density01 gp_el gp_em {
cap drop std_`var'
reghdfe `var', absorb(i.dow##i.admidate_mont##i.finyear i.hyid, savefe) keepsingleton
predict temp_`var', residual
egen std_`var' = std(temp_`var')
}


regress std_adj_docfe std_lsoa_imdscore10 std_lsoa_wa_benefits std_prevyear_admit_em std_prevyear_admit_el std_prevyear_los std_lsoa_density01 std_gp_el std_gp_em 
foreach var in lsoa_imdscore10 lsoa_wa_benefits prevyear_admit_em  prevyear_admit_el prevyear_los lsoa_density01 gp_el gp_em {
di _b[std_`var'] 
matrix observe_`var' = _b[std_`var'] 	
}

** now bootstrap

restore

capture program drop myboot

program define myboot, rclass

* draw my bootstrap sample
preserve
bsample, strata(doctor_id)

cap drop *docfe* *pred_survive30*

reghdfe survive365 c.prevyear_cost, absorb(i.sex##i.derv_age i.black i.mixed i.chinese i.asian i.race_miss i.ynch* i.prevyear_stroke i.di1 i.di2 i.di3 i.di4 i.di5 i.shock i.arythmia i.arthero i.arrest i.dow##i.admidate_mont##i.finyear i.hyid, savefe) keepsingleton
predict residual, residuals

xtreg residual c.std_ami3 c.std_nonami3, fe i(doctor_id)
predict docfe, u
gen sigma_u = e(sigma_u)
gen sigma_e = e(sigma_e)

gen signal_2step = ((sigma_u^2) / ((sigma_u^2) + ((sigma_e^2)/vol)))
gen adj_docfe = docfe*signal_2step

egen std_adj_docfe = std(adj_docfe)

foreach var in lsoa_imdscore10 lsoa_wa_benefits prevyear_admit_em prevyear_admit_el prevyear_los lsoa_density01 gp_el gp_em {
cap drop std_`var'
reghdfe `var', absorb(i.dow##i.admidate_mont##i.finyear i.hyid, savefe) keepsingleton
predict temp_`var', residual
egen std_`var' = std(temp_`var')
}


regress std_adj_docfe std_lsoa_imdscore10 std_lsoa_wa_benefits std_prevyear_admit_em std_prevyear_admit_el std_prevyear_los std_lsoa_density01 std_gp_el std_gp_em
foreach var in lsoa_imdscore10 lsoa_wa_benefits prevyear_admit_em  prevyear_admit_el prevyear_los lsoa_density01 gp_el gp_em{
di _b[std_`var'] 	
return scalar d_`var' = _b[std_`var'] 	
}

restore

end


** Then run bootstrap here

#delimit ;
* Simulate 199 times to test;

simulate diff_lsoa_imdscore10 = r(d_lsoa_imdscore10)
diff_lsoa_wa_benefits = r(d_lsoa_wa_benefits)
diff_prevyear_admit_em = r(d_prevyear_admit_em)
diff_prevyear_admit_el = r(d_prevyear_admit_el)
diff_prevyear_los = r(d_prevyear_los)
diff_lsoa_density01 = r(d_lsoa_density01)
diff_gp_el = r(d_gp_el)
diff_gp_em = r(d_gp_em)
, reps(199) seed(32786105): myboot;

* Boostrap;
bstat, stat(observe_lsoa_imdscore10, observe_lsoa_wa_benefits, observe_prevyear_admit_em, observe_prevyear_admit_el, observe_prevyear_los, observe_lsoa_density01, observe_gp_el, observe_gp_em);


#delimit cr
* Output results
putexcel set "$results/table4_col1.xlsx", replace
putexcel A1 = "Variable"
putexcel B1 = "Coefficient"
putexcel C1 = "Std. error"

matrix b = e(b)'
putexcel A2 = matrix(b), rownames

matrix se = e(se)'
putexcel C2 = matrix(se)


** Then do this for column two (09 onwards sample with extra ED and OP variables)
***********************************************************************************


use "$savedata/masterdata.dta", clear

merge m:1 extract using "$savedata/past_op0517.dta"
keep if _merge==3
drop _merge

merge m:1 extract using "$savedata/pastae0717.dta"
keep if _merge==3
drop _merge

merge m:1 gpprac finyear using "$savedata/gp_pastadmits0517.dta"
drop if _merge==2
drop _merge

* Need to create "last year" variables for use

foreach var in num_cardio_op num_visit_op any_visit_op any_cardio_op{
	
gen prevyear_`var' = `var'2005 if finyear==2006

forval year=2007(1)2017{
	
	local past = `year'-1
	replace prevyear_`var' = `var'`past' if finyear==`year'
}

}


foreach var in any_amb any_cardio_ae any_visit_ae num_amb_ae num_cardio_ae num_visit_ae{
	
gen prevyear_`var' = `var'2007 if finyear==2008

forval year=2009(1)2017{
	
	local past = `year'-1
	replace prevyear_`var' = `var'`past' if finyear==`year'
}

}

drop any_amb* any_visit* any_cardio* num_amb* num_visit* num_cardio*

keep if sample25==1
gen vol = vol25

merge m:1 lsoa using "$inputs/lsoa_areatypes.dta"
drop if _merge==2
drop _merge

* Keep only the data from 2008 onwards as this has ED and OP variables for all
keep if finyear>2008

preserve

reghdfe survive365 c.prevyear_cost, absorb(i.sex##i.derv_age i.black i.mixed i.chinese i.asian i.race_miss i.ynch* i.prevyear_stroke i.di1 i.di2 i.di3 i.di4 i.di5 i.shock i.arythmia i.arthero i.arrest i.dow##i.admidate_mont##i.finyear i.hyid, savefe) keepsingleton
predict residual, residuals

xtreg residual c.std_ami3 c.std_nonami3, fe i(doctor_id)
predict docfe, u
gen sigma_u = e(sigma_u)
gen sigma_e = e(sigma_e)

gen signal_2step = ((sigma_u^2) / ((sigma_u^2) + ((sigma_e^2)/vol)))
gen adj_docfe = docfe*signal_2step

egen std_adj_docfe = std(adj_docfe)

foreach var in lsoa_imdscore10 lsoa_wa_benefits prevyear_admit_em prevyear_admit_el prevyear_los lsoa_density01 prevyear_num_cardio_op prevyear_num_visit_op prevyear_any_visit_op prevyear_any_cardio_op prevyear_any_amb prevyear_any_cardio_ae prevyear_any_visit_ae prevyear_num_amb_ae prevyear_num_cardio_ae prevyear_num_visit_ae gp_el gp_em {
cap drop std_`var'
reghdfe `var', absorb(i.dow##i.admidate_mont##i.finyear i.hyid, savefe) keepsingleton
predict temp_`var', residual
egen std_`var' = std(temp_`var')
}

regress std_adj_docfe std_lsoa_imdscore10 std_lsoa_wa_benefits std_prevyear_admit_em std_prevyear_admit_el std_prevyear_los std_lsoa_density01  std_prevyear_num_cardio_op std_prevyear_num_visit_op  std_prevyear_num_amb_ae std_prevyear_num_cardio_ae std_prevyear_num_visit_ae std_gp_el std_gp_em 
foreach var in lsoa_imdscore10 lsoa_wa_benefits lsoa_density01 prevyear_los prevyear_admit_em prevyear_admit_el gp_el gp_em prevyear_num_visit_op prevyear_num_cardio_op  prevyear_num_visit_ae prevyear_num_cardio_ae prevyear_num_amb_ae {
di _b[std_`var'] 
matrix fobserve_`var' = _b[std_`var'] 	
}


** now bootstrap

restore

capture program drop myboot

program define myboot, rclass

* draw my bootstrap sample
preserve
bsample, strata(doctor_id)

cap drop *docfe* *pred_survive30*

reghdfe survive365 c.prevyear_cost, absorb(i.sex##i.derv_age i.black i.mixed i.chinese i.asian i.race_miss i.ynch* i.prevyear_stroke i.di1 i.di2 i.di3 i.di4 i.di5 i.shock i.arythmia i.arthero i.arrest i.dow##i.admidate_mont##i.finyear i.hyid, savefe) keepsingleton
predict residual, residuals

xtreg residual c.std_ami3 c.std_nonami3, fe i(doctor_id)
predict docfe, u
gen sigma_u = e(sigma_u)
gen sigma_e = e(sigma_e)

gen signal_2step = ((sigma_u^2) / ((sigma_u^2) + ((sigma_e^2)/vol)))
gen adj_docfe = docfe*signal_2step

egen std_adj_docfe = std(adj_docfe)

foreach var in lsoa_imdscore10 lsoa_wa_benefits prevyear_admit_em prevyear_admit_el prevyear_los lsoa_density01 prevyear_num_cardio_op prevyear_num_visit_op prevyear_any_visit_op prevyear_any_cardio_op prevyear_any_amb prevyear_any_cardio_ae prevyear_any_visit_ae prevyear_num_amb_ae prevyear_num_cardio_ae prevyear_num_visit_ae gp_el gp_em {
cap drop std_`var'
reghdfe `var', absorb(i.dow##i.admidate_mont##i.finyear i.hyid, savefe) keepsingleton
predict temp_`var', residual
egen std_`var' = std(temp_`var')
}

regress std_adj_docfe std_lsoa_imdscore10 std_lsoa_wa_benefits std_prevyear_admit_em std_prevyear_admit_el std_prevyear_los std_lsoa_density01 std_prevyear_num_cardio_op std_prevyear_num_visit_op  std_prevyear_num_amb_ae std_prevyear_num_cardio_ae std_prevyear_num_visit_ae std_gp_el std_gp_em 
foreach var in lsoa_imdscore10 lsoa_wa_benefits lsoa_density01 prevyear_los prevyear_admit_em prevyear_admit_el  gp_el gp_em prevyear_num_visit_op prevyear_num_cardio_op  prevyear_num_visit_ae prevyear_num_cardio_ae prevyear_num_amb_ae {
di _b[std_`var'] 
return scalar fd_`var' = _b[std_`var'] 
}

restore

end


** Then run bootstrap here

#delimit ;
* Simulate 199 times to test;

simulate fdiff_lsoa_imdscore10 = r(fd_lsoa_imdscore10)
fdiff_lsoa_wa_benefits = r(fd_lsoa_wa_benefits)
fdiff_lsoa_density01 = r(fd_lsoa_density01)
fdiff_prevyear_los = r(fd_prevyear_los)
fdiff_prevyear_admit_em = r(fd_prevyear_admit_em)
fdiff_prevyear_admit_el = r(fd_prevyear_admit_el)
fdiff_gp_el = r(fd_gp_el)
fdiff_gp_em = r(fd_gp_em)

fdiff_prevyear_num_visit_op  = r(fd_prevyear_num_visit_op)
fdiff_prevyear_num_cardio_op  = r(fd_prevyear_num_cardio_op)
fdiff_prevyear_num_visit_ae  = r(fd_prevyear_num_visit_ae)
fdiff_prevyear_num_cardio_ae  = r(fd_prevyear_num_cardio_ae)
fdiff_prevyear_num_amb_ae  = r(fd_prevyear_num_amb_ae)

, reps(199) seed(32786105): myboot;




#delimit ;
* Boostrap;
bstat, stat(fobserve_lsoa_imdscore10, fobserve_lsoa_wa_benefits, fobserve_lsoa_density01, fobserve_prevyear_los, fobserve_prevyear_admit_em, fobserve_prevyear_admit_el, fobserve_gp_el, fobserve_gp_em,
fobserve_prevyear_num_visit_op, fobserve_prevyear_num_cardio_op, fobserve_prevyear_num_visit_ae, fobserve_prevyear_num_cardio_ae, fobserve_prevyear_num_amb_ae 
);


#delimit cr
* Output results
putexcel set "$results/table4_col2.xlsx", replace
putexcel A1 = "Variable"
putexcel B1 = "Coefficient"
putexcel C1 = "Std. error"

matrix b = e(b)'
putexcel A2 = matrix(b), rownames

matrix se = e(se)'
putexcel C2 = matrix(se)





